module fft 
use technical
use omplib
implicit none
include 'fftw_f77.i'
include 'fft.h'
integer*8,save :: plan1d,iplan1d,plan2d,iplan2d
integer,save :: nfft1,nfft2,nthreads
contains
!======================
  subroutine initialize_fft_1dcomplex2real(N)
  integer :: N
!
  call rfftw_f77_create_plan(iplan1d,N,FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE)
!
  endsubroutine initialize_fft_1dcomplex2real 
!======================
  subroutine initialize_fft_1dreal2complex(N)
  integer :: N
!
  call rfftw_f77_create_plan(plan1d,N,FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE)
!
  endsubroutine initialize_fft_1dreal2complex
!======================
  subroutine initialize_2dfft(dim,lfft_in_place)
  integer,dimension(2) :: dim
  logical :: lfft_in_place
  integer :: success 

!fftw/fftwnd/rfftw/rfftwnd replaced by fftw_f77/fftwnd_f77/rfftw_f77/rfftwnd_f77.
  
  call fftw_f77_threads_init(success) 
  if (success.eq.0) then
    nfft1=dim(1)
    nfft2=dim(2)
    call initialize_fft_2dcomplex2real(dim,lfft_in_place)
    call initialize_fft_2dreal2complex(dim,lfft_in_place)
!$OMP PARALLEL 
!$OMP MASTER
    call get_total_threads(nthreads) 
    write(*,*) 'Running FFTW with threads, no. of threads=',nthreads
!$OMP END MASTER
!$OMP END PARALLEL 
  else
    call fatal_error('Initializing FFTW','Could not initialize threaded FFTW')
  endif 
!
   endsubroutine  initialize_2dfft
!======================
  subroutine initialize_fft_2dcomplex2real(dim,lfft_in_place)
  integer,dimension(2) :: dim
  logical :: lfft_in_place
!
  if (lfft_in_place) then
    call rfftwnd_f77_create_plan(iplan2d,2,dim,FFTW_COMPLEX_TO_REAL, &
         FFTW_MEASURE + FFTW_IN_PLACE)
  else
    call rfftwnd_f77_create_plan(iplan2d,2,dim,FFTW_COMPLEX_TO_REAL, &
         FFTW_MEASURE)
  endif
!
  endsubroutine initialize_fft_2dcomplex2real 
!======================
  subroutine initialize_fft_2dreal2complex(dim,lfft_in_place)
  integer,dimension(2) :: dim
  logical :: lfft_in_place
!
  if (lfft_in_place) then
    call rfftwnd_f77_create_plan(plan2d,2,dim,FFTW_REAL_TO_COMPLEX, &
         FFTW_MEASURE + FFTW_IN_PLACE)
  else
    call rfftwnd_f77_create_plan(plan2d,2,dim,FFTW_REAL_TO_COMPLEX, &
         FFTW_MEASURE)
  endif
!
  endsubroutine initialize_fft_2dreal2complex
!======================
  subroutine fft1d_complex2real(N,gk,gx)
  integer :: N
  double precision,dimension(N) :: gk,gx
!
  call rfftw_f77_one(iplan1d,gk,gx)
!
  endsubroutine fft1d_complex2real
!======================
  subroutine fft1d_real2complex(N,gx,gk)
  integer :: N
  double precision,dimension(N) :: gx,gk
!
  call rfftw_f77_one(plan1d,gx,gk)
!
  endsubroutine fft1d_real2complex
!============================================
  subroutine fft_2d_complex2real(gk)
  double precision,dimension(nfft1+2,nfft2) :: gk
!
    call rfftwnd_f77_threads_one_complex_to_real(nthreads,iplan2d,gk,0)

  endsubroutine fft_2d_complex2real
!============================================
  subroutine fft_2d_real2complex(gx)
  double precision,dimension(nfft1+2,nfft2) :: gx
!
    call rfftwnd_f77_threads_one_real_to_complex(nthreads,plan2d,gx,0)

  endsubroutine fft_2d_real2complex
!============================================
  subroutine fft_2d_cleanup()
  call rfftwnd_f77_destroy_plan(iplan2d)
  call rfftwnd_f77_destroy_plan(plan2d)
  endsubroutine fft_2d_cleanup
!============================================
end module fft
